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Abstract 

We have studied the superconducting-insulating phase transition in a disordered two-dimensional 
Josephson junction array, using quantum Monte Carlo techniques. We consider disorder in both 
the capacitive energies and in the values of the offset charges. The calculated phase diagram shows 
that the lobe structure of the phase diagram disappears for sufficiently strong disorder in the offset 
charge. Our results agree quite well with previous calculations carried out using a mean-field 
approximation . 



PACS numbers: 74.25.Dw, 05.30.Jp, 85.25.Cp 



I. INTRODUCTION 



A Josephson junction array (JJA) consists of a collection of superconducting islands 
connected by Josephson coupling. The coupling can arise from tunnel junctions through 
an insulating layer, or from the proximity effect. The arrays themselves can be produced 
experimentally in a wide range of geometries, and with a great variety of individual junc- 
tion parameters. They can serve as valuable model systems for studying quantum phase 
transitions^ under conditions such that the experimental parameters can be readily tuned 3 . 
Recently, it has also been proposed that small groups of Josephson junctions may serve as 
quantum bits (qubits) in quantum information technology. In this case, the quantum logic 
operations are performed by manipulating experimental parameters such as gate voltages or 
magnetic fields^. 

If the superconducting islands in a Josephson array are sufficiently large, the array is 
believed to become superconducting in two stages. First, at a temperature near the bulk 
transition temperature, each island becomes superconducting. Secondly, at a lower temper- 
ature, and provided that the array is at least two-dimensional, thermal fluctuations become 
sufficiently weak that a global phase coherence can be established throughout the array^. 

If, however, the islands in the array are small, such phase coherence may not be estab- 
lished even at temperature T = 0. The reason is that there is a second energy which becomes 
important in small grains, namely the charging energy of the grains^. The crucial physics 
is then determined by the competition between this charging energy and the Josephson 
coupling energy. If the charging energy is sufficiently large, then it becomes prohibitively 
expensive energetically to transfer Cooper pairs from grain to grain, and the array becomes 
an insulator, even though each grain is in its superconducting state. If, on the other hand, 
the Josephson dominates, the array becomes phase coherent at T = 0. In this case, Cooper 
pairs can tunnel between neighboring grains and the array as a whole will be in the super- 
conducting phase. Thus, by tuning the ratio of the charging and Josephson energies, one 
can cause the array to undergo a quantum phase transition between a superconducting and 
an insulating stated 3 -. 

In the limit of very large number of Cooper pairs per grain, the Josephson junction array 
with diagonal charging energy is equivalent to the Bose Hubbard model (BHM), which 
describes soft core bosons hopping on a lattice with on-site Coulomb interactions. The 
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BHM has previously been extensively studied by Fisher etal — , who constructed its T = 
phase diagram. In the absence of disorder, the phase diagram shows two phases: a Mott 
insulating phase and a superfluid phase. In the presence of disorder, an additional phase, 
known as the Bose glass phase, emerges^. Recently, the phase diagram of the disordered 
BHM in 2-dimensions (2D) has been studied by Lee et al using QMC3. 

The aim of this paper is to construct the phase diagram for the superconducting-insulating 
(SI) transition in a particular model for a JJA, including the effects of disorder. Such 
disorder is clearly unavoidable in most practical systems. The model Hamiltonian we study 
has previously been investigated by several authors in ordered and disordered arrays, using 
both mean-field theory (MFT) and quantum Monte Carlo (QMC) techniques. A central goal 
of our work is to compare the MFT of Refs.^ and QMC results in a disordered array in order 
to check the accuracy of the MFT. We shall show that the MFT is generally satisfactory, 
even for disordered systems. 

The remainder of this paper is organized as follows. In the next section, we present the 
model Hamiltonian and describe the QMC algorithm we use. Our numerical results are 
presented in Section III, followed by a brief concluding discussion in Section IV. 

II. MODEL HAMILTONIAN AND QMC ALGORITHM 

Our model of a JJA Hamiltonian involves two types of degrees of freedom: the number 
of excess Cooper pairs hi on the ith grain, and the phase fa of the superconducting order 
parameter on the ith grain, hi and (pi are taken to be quantum-mechanically conjugate 
variables with commutation relations [hi, fa] = —iSij. We consider the following model 
Hamiltonian on a square lattice in 2D: 

H = \ U ii(^i ~ ^if ~ E J £ c Os(0i - fa). (1) 

Here Ej is the Josephson coupling strength between nearest neighbors denoted by 
(assumed to be the same for all nearest neighbor pairs), and Uu is the charging energy of 
the ith grain. We expect that Uu = q 2 /Cu, where Cu represents the capacitance of the ith 
grain with respect to ground, and q = 2e is the charge of a Cooper pair. In a more general 
model, the charging part of the energy would be written ^Uij(hi—hi)(hj—hj), where U^ is an 
element of the charging energy matrix. In earlier calculations using MFT, nearest-neighbor 
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and next-nearest-neighbor terms have been included, but these are numerically more difficult 
to include in QMC than the diagonal terms. Finally, the quantity fii represents the "offset 
charge." is related to the voltage between the zth grain and a common ground plane. 
In an ordered array, fii should be independent of i. We omit any dissipative terms arising 
from Ohmic shunts in the junctions, and we assume that amplitude fluctuations of the gap 
on the individual grains can be neglected; this neglect of amplitude fluctuations should be 
reasonable if T -C T^, the single-grain transition temperature. Thus tunneling of charge 
between grains involves only Cooper pairs, not single electrons. 

In our calculations, we include disorder in two terms in the Hamiltonian: the diagonal 
charging element Ua and the offset charge n,. Disorder in Ua may arise from randomness in 
the size of the individual grains, while disorder in fii could arise from random offset charges 
near the superconducting grains. Such disorder is unavoidable in practical realizations of 
Josephson arrays. 

It is worthwhile to discuss the effects of the offset charges fii qualitatively. Such charges 
can arise in two ways: (i) from a voltage applied between the array and the substrate, and 
(ii) from random charges in the substrate. In the latter case, n« is a random variable. In 
arrays with only a few grains, fii can, in principle, be tuned individually to desired values. 
However, such tuning, although possibly still achievable in principle, may be difficult to 
attain in practice in arrays with many junctions. 

In ordered arrays, the most important effect of the offset charge is to reduce the region of 
the phase diagram where the Mott insulating phase is stable. The reason is that, if n« 7^ 0, 
charge fluctuations become energetically less expensive; it is therefore correspondingly easier 
to establish phase coherence and to destroy the Mott insulator. 

For the disordered array, we determine the phase diagram of the Hamiltonian (0) using 
QMC techniquesii>^, and compare the results with MFT 7 . Following the method of Ref.— , 
we first map the 2D Hamiltonian Hamiltonian 7i [Eq. (JI|) ]. onto a classical 3D Hamiltonian. 
This is accomplished by the standard procedure of transforming the partition function Z = 
Tre~P' H , where (3 = l/ksT, into an Euclidean path integral along the imaginary time axis 
r from to Tif3. To do the transformation, one breaks up the time integral into L T small 
time steps, each of length e = h/3/L T , and uses the identity e~^ n = e~ eH . . .e~ eH . Next, a 
complete set of eigenstates (for example, the eigenstates of the Cooper pair number operator 
for the zth grain) are inserted at each imaginary time step. The next step is to rewrite the 
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Josephson term using the Villain approximation^ as e~ xcos ^ ~ J2m=-oc e _a: ^^^ _27rm ^ / 2 , 
where f(x) = {2x\ti[Iq(x)/Ii(x)}}~ 1 , and J„(x) is the modified Bessel function of order n. 
For large values of x, f(x) ~ 1. 

Finally, the partition function Z can be expressed as Z = Tre~ s , where the action 5 is 
defined as 

S = o 6 E ^ ( J ft } - ^,r)(^> - %r) 



E 1 4? I 2 - (2) 



2Ejef(Eje, 

In Eq. (J2J), the new degrees of freedom are integer- valued current loops jj^^a — x, y), which 
live on the nodes of the 3-dimensional (xyr) lattice, and satisfy the continuity equation 
J2 a =x,y,d Q Ji°r > + 9 T J^J = at every lattice point. The time components of the current 
operators, jf^, represent the Cooper pair number operators n itT along their world lines. 
Note that S is a classical action in (d + 1) dimensions while the original Ti is a quantum 
Hamiltonian in d dimensions. 

To evaluate S, we use the standard Metropolis algorithm to generate configurations of 
the currents Jj )T at inverse temperature (3. It is convenient to work in the grand canonical 
ensemble. The system size is assumed to be a parallelepiped with dimensions (L, L, L T ) 
along the space and imaginary time axes respectively, with periodic boundary conditions in 
all three directions. 

We locate the superconducting-insulating (SI) phase transition by examining the super- 
fluid density p, which is proportional to the stiffness of the system against a twist in the 
phase. As has been shown previously^, p can be expressed in terms of the current variables 

p = jk; {{w{x))) - (3) 



Here 



J2i, T Ji^- ^ S ^ ne sc, - ca ll e d winding number in the x-direction, and ((•••)) denotes 
both a grand canonical and a disorder average. If p is non-zero in the thermodynamic limit 
of a large system, then there is long-range phase coherence and the system will be in the 
superconducting phase. 
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FIG. 1: Phase diagram for a model Josephson junction array described by the Hamiltonian (1), 
plotted as a function of the offset charge ft as calculated using QMC and MF theory, in the absence 
of disorder. The curves show the critical value of J/Uq separating the superconducting (S) phase 
from the insulating (I) phase. The temperature is ksT = 0.03U and the QMC calculations are 
done for a lattice of size 20 3 . The line segments simply connect the calculated points 

III. RESULTS 

We have carried out our simulations for three different system sizes: 8 3 , 10 3 , and 14 3 , 
and have usually averaged over ~ 100 — 500 realizations of the disorder. For system with 
no disorder we studied system with a size of 20 3 . For each disorder realization, the system 
is equilibrated by slowly annealing in temperature, over about 2000 x L 2 L T passes through 
the entire system for each value of Ej considered, followed by approximately twice as many 
steps over which the grand canonical averages are computed. 

Figure n shows the phase diagram for a JJA at k^T = 0.3U with diagonal charging 
energy and no disorder, as calculated using QMC as described above. For reference, we also 
show the same phase diagram as calculated using MFT&£. The regions of the phase diagram 
denoted S or / correspond to the superconducting and Mott insulating regions, respectively. 
The two curves show the phase boundaries as calculated using MFT (squares) and QMC 
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FIG. 2: (a) Phase diagram for the model Hamiltonian (1) as calculated using QMC (upper panel) 
and MFT theory (lower panel), including the effects of disorder. The offset charge is chosen from a 
uniform distribution between n — a and n + a. For the QMC calculation, we use a 14 3 lattice and 
fix the temperature at ksT = 0.03[/o. (b) Same as (a) except that the offset charge is chosen from 
a Gaussian distribution of standard deviation a centered at n. In these calculations, line segments 
connect calculated points. 

(triangles). Although in this and later Figures we show results for < n < 1, calculations 
were actually carried out only in the region < n < 0.5, since the phase diagram must be 
symmetric about n = 0.5. Our findings here are similar to what is found by Otterlo etal in 
RefA 

We have considered two types of disorder: randomness in the offset charges n,, and 
randomness in the diagonal charging energies Uu. Disorder in fii is the analog of chemical 
potential disorder in the BHM, while disorder in the Ua's is analogous to randomness in the 
mean-field on-site Coulomb energy in the BHM. 

We begin by describing our results for offset charge randomness. We have carried out 
calculations with two different random distributions of n's. In the first case, we choose the 
offset charge at each lattice site so that fii = n + <7i, where (7j is a random number uniformly 
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FIG. 3: Comparison of the phase diagram as calculated using QMC (triangles) and MF the- 
ory (squares), assuming a uniform distribution of n and several values of the disorder strength 
parameter a, as indicated in each plot. Other parameters are the same as in Fig. [21 

distributed in the range [— a, a] , where o represents the strength of disorder. In the second 
case, we choose fii to have a Gaussian distribution with a mean value of n and standard 
deviation a, i.e., with a probability distribution proportional to e -( n i- n ) 2 / 2 °" 2 . i n both cases, 
the values of n« on different sites are taken as uncorrelated. 

In Fig. |21 (a), we show the phase diagram for a disordered JJA, assuming a flat distribution 
of Hi, and as calculated using QMC (upper panel), and MFT (lower panel). In each case, 
the different curves correspond to the phase diagrams for different disorder strengths, as 
indicated in each Figure. In Fig. 121(b), we show the corresponding results, but for Gaussian 
disorder with various standard deviations a, as indicated in the figure. 

Several features are noticeable from the QMC calculations. First, at values of n close to 
or 1, disorder reduces the stability region of the insulating phase, whereas near n = 1/2, 
disorder increases the region in which the insulating phase is stable. These effects are 
readily understood: a finite disorder means that one is effectively calculating the phase 
diagram over some average region of n, and therefore the sharp lobes seen in the phase 
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FIG. 4: Same as Fig. |21 except that the offset charge n is chosen from a Gaussian distribution of 
standard deviation a. 

diagram for non-random n should become less distinct as the disorder increases. Indeed, 
the lobe structure shrinks and tends to disappear for large values of disorder. A similar 
effect was also emphasized by Fisher etal for the Bose Hubbard modeR This behavior is 
not surprising, because in the limit of large disorder strength a, the Hamiltonian becomes 
independent of n: all values of n are therefore equivalent. (This equivalence is exact for the 
uniform distribution at o = 0.5.) 

Results from the QMC and MFT calculations are compared in a different way in Figs. E3 
and H]for the two kinds of disorder in nj. As is evident from both Figures, the MFT results 
agree very well with those from QMC results for both types of disorder, except possibly near 
the tips of the lobes. Near n = 0.5, the critical value of J/Uq is slightly underestimated by 
MFT. A similar underestimate occurs for large values of the disorder strength a. In this 
case, MFT underestimates the critical value of J /Uq by about 10%. 

Figure El shows the results of calculations in which disorder is included in the diagonal 
charging energy but not in the n^'s. Specifically, we choose Uu at random from values 
uniformly distributed in the interval [Uo(l — A),t/ (1 + A)]. We have considered several 
values of the disorder strength parameter A, as indicated in the legend. Clearly, randomness 
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FIG. 5: Phase diagram for a disordered 2D Josephson junction array described by the Hamiltonian 
(1), with charging energy disorder. Only the diagonal elements Uu of the charging energy are 
nonzero; the Ua's are distributed randomly and uniformly in the interval [Uq(1 — A), Uq(1 + A)]. 
The values of the disorder parameter A are shown in the Figure. Other parameters are the same 
as in Fig. El 

in the Uu always increases the region of the phase diagram in which the S phase is stable, 
whatever the value of n. This trend is not surprising, because fluctuations in the Ua's 
lead to greater fluctuations in the number of Cooper pairs on each island. Thus, by the 
uncertainty principle, the phase fluctuations should be relatively reduced, and therefore the 
S (phase-ordered) state should become more favorable, as seen in our calculations. 

IV. DISCUSSION AND SUMMARY 

In this paper, we have studied the phase diagram of a model Josephson junction array 
at low temperatures using quantum Monte Carlo techniques. Our goal was to analyze the 
effects of disorder in both the offset charges fii and grain sizes (which affect the diagonal 
elements of the charging energy matrix, Uu). For disorder in the n^s, we find that disorder 
favors the S phase when the average offset charge n is close to an integer number of Cooper 
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pairs, but that it favors the I phase near n = 1/2. We also find that for large disorder, 
the average value n has little effect on the phase diagram at all. As for the second kind 
of disorder, we find that disorder in the L^j's tends to favor the S phase for all values of 
n, as shown in Fig. 5. Finally, we have found that the critical value of Ej/U at the phase 
boundary, as calculated from the MFT of Grignani et al — , agrees well with that calculated 
using QMC (almost perfectly for small disorder and within 10% even for large disorder). 
Thus, the MFT gives good results for most of the phase diagram. 
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